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In this third and final paper of a series, elastic properties of numerically simulated isotropic 
packings of spherical beads assembled by different procedures, as described in the first companion 
paper, and then subjected to a varying confining pressure, as reported in the second companion 
paper, are investigated. In addition to the pressure, which determines the stiffness of contacts 
because of Hertz's law, elastic moduli are chiefly sensitive to the coordination number z, which should 
not be regarded as a function of the packing density. Comparisons of numerical and experimental 
results for glass beads in the lOkPa — lOMPa pressure range reveal similar differences between 
dry samples prepared in a dense state by vibrations and lubricated packings, so that the greater 
stiffness of the latter, in spite of their lower density, can be attributed to a larger coordination 
number. Effective medium type approaches, or Voigt and Reuss bounds, provide good estimates of 
bulk modulus B, which can be accurately bracketed, but badly fail for shear modulus G, especially 
in low 2 configurations under low pressure. This is due to the different response of tenuous, fragile 
networks to changes in load direction, as compared to load intensity. In poorly coordinated packings, 
the shear modulus, normalized by the average contact stiffness, tends to vary proportionally to the 
degree of force indeterminacy per unit volume, even though this quantity does not vanish in the 
rigid limit. The elastic range extends to small strain intervals and compares well with experimental 
observations on sands. The origins of nonelastic response are discussed. We conclude that elastic 
moduli provide access to mechanically important information about coordination numbers, which 
escape direct measurement techniques, and indicate further perspectives. 

PACS numbers: 45.70.-n, 83.80.Fg, 46.65.+g, 62.20. Fe 



I. INTRODUCTION 

The mechanical properties of granular materials and 
their relations to the packing microstructure are cur- 
rently being investigated by many research groups. As a 
simple model, long studied for its geometric aspects [1,|31j 
the packing of equal-sized spherical balls is also mechan- 
ically characterized in the laboratory [1, 0, 0, B j and 
by numerical means, relying on discrete, granular level 
modehng [1, [M [O, 111, 111, 13 ■ 

The present paper is the last one in a series of three, 
about geometric and mechanical properties of bead pack- 
ings obtained by numerical simulations. It focusses on 
the elastic properties of isotropically compressed sam- 
ples. The study is based on the configurations for which 
the packing processes and resulting microstructure were 
studied in paper I while paper II 2] reported on the 
effects of isotropic compressions and pressure cycles. The 
results presented here, although perhaps better appreci- 
ated on knowing the contents of papers I and II, can be 
understood without reading these previous contributions. 
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Elastic properties of granular assemblies are probed 
when small stress increments are superimposed on a pre- 
strcsscd equilibrium configuration, either on controlling 
very small strains in a static experiment [isl [l^.[l7l or 
in dynamical ones, relying on resonance modes llSl . Il9l . 
[13 , or sound propagation 0, 0, [H [13, H HH US, IM IM ■ 
Elastic behavior of granular materials is only applica- 
ble for very small strain increments, typically of or- 
der 10~^ or even 10^^ in usual conditions, i.e., with 
sands under confining stresses between 10 kPa and a few 
MPa [6, 15, 16, 17]. It has been checked in such cases that 
static measurements of elastic moduli, with devices accu- 
rate enough to control such small strains are consistent 
with "dynamical" ones, i.e. deduced from experiments 
on wave propagation or resonance frequencies. Experi- 
mental soil mechanics have achieved a high level of so- 
phistication, with significant progress over the last twenty 
years (25l . [2a | , and accurate measurements of the mechan- 
ical response of granular materials in the very small strain 
regime are one example thereof. Coincidence of elastic 
moduli values obtained by different means is reported, 
e.g., in [H, [13,113. 

The elastic moduli should not be confused with the 
slopes of stress-strain curves on the scale of the strain 
level (usually in the 1% range) corresponding to the full 
mobilization of internal friction. Such slopes are consid- 
erably smaller than true elastic moduli (by more than an 
order of magnitude), and do not correspond to an clas- 
tic, reversible response. In this respect, the frequent use, 
for engineering applications, of a simplified elastoplastic 
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behavior, in which the material is linear elastic until the 
Mohr-Coulomb criterion for plasticity is reached, as pre- 
sented in Ref. [l^], should not be misinterpreted. Such 
crude models, in which strains are elastic and reversible 
up to the complete mobilization of internal friction, are 
resorted to in engineering practice when detailed infor- 
mation on the constitutive law are not available, but the 
"elastic moduli" introduced in those simplified constitu- 
tive relations are merely convenient parameters enabling 
one to perform approximate calculations. 

In micromcchanical [28, 29] and numerical [1, [13, [13] 
studies, elastic properties are associated with the defor- 
mations of a fixed contact network, and should therefore 
correspond to the "true elastic" behavior observed in the 
laboratory for very small strain intervals. Indeed, ex- 
cept in very special situations in which the effects of fric- 
tion are suppressed and geometric rcstructuration is re- 
versible 131, ^52] , the irreversible changes associated with 
network alterations or rearrangements preclude all kind 
of elastic modelling. Elastic properties are attached to 
one specific contact set, and hence of limited relevance to 
the rheology of granular materials. Nervertheless, elas- 
tic properties are interesting because they might provide 
access, in a non-destructive way, to geometric data on 
the contact network, such as coordination numbers. Such 
variables are still virtually inaccessible to direct measure- 
ments, even with sophisticated visualization techniques, 
as emphasized in paper I but they are very likely, in 
turn, to influence the constitutive laws for larger strains. 

This paper is organized in the following way. We first 
recall the properties of the model material we are study- 
ing (Section |Tl| , along with basic definitions and prop- 
erties pertaining to the elasticity of granular packings. 
Then, useful results on the pressure-dependent internal 
states of the various types of configurations introduced 
and studied in papers I and II [ll, I2l are summarized in 
Section UTTl Next, the values of elastic constants in the 
different configuration series, as a function of (isotropic) 
confining pressure, are presented in Section IIVI where 
their relations to internal structure are also discussed. 
Section |V] is devoted to the particular behavior of elastic 
moduli in the tenuous contact networks of poorly coor- 
dinated configurations. Numerical results are confronted 
to experimental ones in Section IVII Some results about 
the extension of the elastic range are given in Section lVlII 
Section IVlIII discusses the results and indicates some fur- 
ther perspectives. 



we request all three diagonal components of the Cauchy 
stress tensor, denoted as Gaa to be equal, in equilibrium, 
to prescribed values S^, all chosen to coincide with a 
pressure P in this study of isotropic states. Differences 
between aaa and Sq, entail some evolution in the cell size 
parameters. Uaa is given in equilibrium by the classical 
formula: 



i<3 



(1) 



where the sum runs over all pairs in contact, F^"' is the 
a coordinate of the force exerted by grain i onto its 

neighbor j at their contact, while r^^"^ is the a coordinate 
of vector Vij^ pointing from the center of i to the center 
of j. From Eqn. ([T]) one can easily deduce a simple and 
useful relation between pressure P = (cth -I- 022 + o'33)/3 
and the average normal contact force {N) between mono- 
sized spheres of diameter a, involving solid fraction $ and 
coordination number z : 



P 



z$(iV) 



(2) 



The corresponding dynamical equations used to impose 
stresses are described in paper I in, and we briefiy recall 
here the essential ingredients of the model for a study 
of the elastic response of packings that have been first 
assembled and compressed, as reported in papers I and 
II. Dynamical aspects of the model, in particular (inertia, 
viscous dissipation) play no role in the determination of 
elastic moduli, for which our calculations are based on the 
building of the stiffness matrix of the contact network. 



A. Local stiffnesses 

We consider spherical beads of diameter a, with the 
elastic properties of glass: Young modulus E — 70 GPa, 
Poisson ratio ly = 0.3. They interact in their contacts by 
the Hertz law, which states that the elastic normal force 
N is proportional to h^^'^, h being the normal deflection 
of the contact, so that the incremental normal stiffness 
with the notation E = jS^, is given by: 

dh 2 2 ^ ' 



II. NUMERICAL MODEL AND BASIC 
DEFINITIONS 

Packings of n spherical beads are simulated with 
molecular dynamics, in which equations of motion re- 
sulting from Newton's laws are solved for the particle 
positions and rotations. Thanks to a suitably adapted 
form of the Parrinello- Rahman deformable cell molecular 
dynamics technique [H, [13] , as described in paper I [l| , 



The tangential elastic force is to be incrementally evalu- 
ated with a simplified Mindlin-Deresiewicz form [ssj on 
assuming the tangential stiffness Kt to stay proportional 
to Kn, and hence a function of N: 

Kt = arK^ with ax = . (4) 

2 — 

The tangential force T is constrainted by the Coulomb 
condition ]]T]] < iiN, with friction coefficient ^ set to 
0.3, and additional conditions are introduced to ensure 
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thermodynamic consistency and objectivity (see paper 
I [U, Sec. II, appendices A and B]). 

All our results will be stated in a form independent of 
bead diameter a. All dimensionless results - such as e.g. 
ratios of elastic moduli of the granular material to E - 
depend, in addition to on a reduced stiffness parameter 
we define as 

2/3 



Typical values of contact deflections h scale as K~^a. 

Let Nc denote the number of force-carrying contacts. 
In every contacting pair i-j, we arbitrarily choose a 
"first" grain i and a "second" one j, and define the rel- 
ative dispacement vector Ju^ as the difference between 
the displacement of the contact point as belonging to 
solid i and its displacement as a point belonging to solid 
j, both regarded as rigid. 

In each contact the force F.^ that is transmitted from 
i to j is split into its normal and tangential components 
as Fij — NijXiij + Tij. The static contact law relates 
the SA^c-dimensional contact force increment vector Af , 
formed with the values ANij, AT^ of the normal and 
tangential parts of all contact force increments, to Su: 



are degrees of freedom of the system, while all three di- 
agonal components of the stress tensor are externally im- 
posed We use three strain parameters defined as 
the relative changes of those lengths, from their values in 
a reference state: 

With our conventions, shrinking strains and compressive 
stresses are positive. 

Let us now recall the definition of the rigidity ma- 
trix (not to be confused with the stiffness matrix), as 
introduced in paper I. The grain center displacements 
(ui)i<i<n are conveniently written as 



(7) 



with a set of displacements Uj satisfying periodic bound- 
ary conditions in the cell with the current dimensions, 
e denoting the diagonal strain matrix with coefficients 
Eq . Gathering all coordinates of particle (periodic) dis- 
placements and rotation increments, and strain parame- 
ters one defines a displacement vector in a space with 
dimension equal to the number of degrees of freedom 
Nf = 3n -I- 3 (recall n is the number of beads). 



U= ((u,,Aa 



l<a<3 



(8) 



Af /C • 5u. 



(5) 



This defines the [iN^. x iN^) matrix of contact stiffnesses 
K. JC is block diagonal (it does not couple different con- 
tacts), and we shall refer to it as the local stiffness matrix. 
The 3x3 block of JC corresponding to contact z, j, JC^^ 
is diagonal itself provided friction is not fully mobilized, 
and contains stiffnesses KN{hij) and (twice in 3 dimen- 
sions) Krihij) as given by Q and (j4]): 

















(6) 



This simple form of IC.^. ignores some non-diagonal terms 
that appear when the normal force decreases, or for slid- 
ing contacts. The effects of those terms are discussed and 
tested in Appendix El as well as the possible use of more 
elaborate contact laws [1^ , in which Kt / Kn is not kept 
constant. 



B. Global rigidity and stiffness matrices 

When elastic properties are investigated, small dis- 
placements about an equilibrium configuration are dealt 
with to first order (as an infinitesimal motion, i.e. just 
like velocities), and related to small increments of applied 
forces, moments and stresses. We use periodic bound- 
ary conditions in our simulations, and the dimensions 
La [a = I, 2, 3) of the parallelipipedic simulation cell 



The normal unit vector riy points from i to j (along 
the line joining centers for spheres). The relative dis- 
placement 6uij , for spherical grains with radius R, reads 



Su, 



X Rui 



59 j X Rxiij 



(9) 



in which r,y is the vector pointing from the center of the 
first sphere i to the nearest image (by the periodic trans- 
lation group of the boundary conditions) of the center of 
the second one j. The normal part <5u^ of (5u.y is the 
increment of normal deflection hij in the contact. 

The rigidity matrix G is 3Nc x iV/-dimensional, it is 
defined by the linear correspondence expressed by rela- 
tion ([9]), which transforms U into the 3iVc-dimensional 
vector of relative displacements at contacts Su: 



Su = G \J 



(10) 



External forces F^ and moments Ti applied to grain 
centers, and diagonal Cauchy stress components can 
be gathered in one A^/-dimensional load vector F"'^*: 



((F„r,) 

1< t<n 1 



(11) 



chosen such that the work in a small motion is equal to 
poxt .-(j^ The equilibrium equations - the statements that 
contact forces f balance load F"''* ~ is simply written with 
the tranposed rigidity matrix, as 



(12) 



This is of course easily checked on writing down all force 
and moment coordinates, as well as the equilibrium form 
of stresses, Eqn. [T]with aaa = Sq, for all a. 
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Given (jlOp . which defines the rigidity matrix, relation 
(fT^ is a statement of the theorem of virtual work (see 
paper I [H). 

Returning to the case of small displacements associated 
with a load increment AF'^^', one may write, to first 
order in U, 

^pext = ^ . U, (13) 

with a total stiffness matrix K, comprising two parts, 
K^^-* and K*-^-*, which we respectively refer to as the con- 
stitutive and geometric stiffness matrices. K^^^ results 
from Eqns . [ini El and [m 

l< = • £ • £ (14) 

K^^-* is due to the change of the geometry of the packing, 
and is written down in AppcndixiBl where it is also shown 
to be negligible in general. 

To the rigidity matrix are associated the concepts of 
force and velocity (or displacement) indeterminacy, of 
relative displacement compatibility and of static admis- 
sibility of contact forces. Those definitions are given in 
paper I, where they are used to discuss the limit of iso- 
static packings [s^l- 

C. Grain-level and macroscopic elasticity 

We now discuss the conditions for which the response 
to load increments of a prestressed granular packing in 
mechanical equilibrium can be described as elastic, and 
explain how macroscopic elastic moduli are computed in 
our simulations. 



1. Some necessary approximations 

Elasticity implies the existence of an elastic potential 
(a function of U) from which forces are derived. If force 
increments are written as linearly depending on displace- 
ments, as in (jl3p . the corresponding stiffness matrix K 
should be unique (the same for all U vectors) and sym- 
metric. Strictly speaking, the differences in the form 
of the local stiffness matrix IC on the direction of dis- 
placements precludes the definition of an elastic response. 
This is discussed in Appendix [A] We could check that 
this effect is quantitatively negligible, and that taking all 
)C . block values as in Eqn. ([6|) is a very good approxima- 

— ^3 

tion. Likewise, the implementation of more complicated 
contact laws accounting for a gradual change of tangen- 
tial stiffness Kt as a function of ||T|| (see Appendix [X}, 
entails negligible changes in the values of elastic moduli. 

Moreover, the non-symmetric geometric contribution 
K^^-* is also negligible - provided simple precautions are 
applied to prevent the free motion of divalent beads, as 
explained in Appendix [B) These motions, as checked in 
paper I do not jeopardize global stability, and they 



are the only mechanisms (i.e., non-zero elements of the 
kernel of rigidity matrix G) of the backbone. On us- 
ing the symmetric diagonal form ([6]) for all contacts, IC 
will be symmetric and positive definite. In view of (jl4[) . 
the kernel of K'"'^' (the "floppy modes" of the constitu- 
tive stiffness matrix) coincides with the kernel of G (the 
"mechanisms"). A suitable elimination of the localized 
free motion of 2-coordinated grains (see Appendix [B|) en- 
ables one to work with a positive definite global stiffness 
matrix, and therefore with a well-behaved elastic poten- 
tial energy. The elastic regime, however, is only defined 
on approximating the contact laws as elastic. 

Finally, as explained in Appendix [Ul initial confining 
stresses entail very small corrections to elastic moduli, of 
order P, which we also neglect. 

2. Computation of elastic moduli 

In order to evaluate macroscopic elastic moduli or com- 
pliances, one can apply stress increments and measure 
the resulting strains. With our choice of boundary con- 
ditions and degrees of freedom, we choose load incre- 
ments AF*''^' with all coordinates set to zero except one 
of the three last ones, say QAaaa, corresponding to a 
diagonal stress increment, according to definition (jlip . 
Then we solve the system of equations for the un- 
known displacement vector U. Its 3 last coordinates 
are identified as diagonal strain components, according 
to definition ([5]). The effective elastic properties of the 
packing being isotropic, we obtain = Ua/E*, and 
ep = —v*UalE* for /9 ^ a, in which E* and v* are the 
effective macroscopic Young modulus and Poisson coef- 
ficient of the bead packing. On changing a one obtains 
different estimates of those macroscopic properties, which 
should coincide in the limit of large systems. 

This procedure is used with the numerical samples pre- 
pared in equilibrium states under varying isotropic pres- 
sure P, as explained in papers I and II (see Section IIIII 
below). All results are averaged over sets of 5 statisti- 
cally similar samples of n = 4000 grains each, and error 
bars on curves (which are often as small as the symbols) 
correspond to one standard deviation on each side of the 
mean value. 



3. Minimization properties 

In order to write down bounds on macroscopic elastic 
moduli, the following minimization properties are use- 
ful [131 ■ First, solving for U is equivalent to mini- 
mizing the following potential energy: 

W'i(U) = iu-K-U- AF'''^' -U. (15) 

Then, the contact force vector increment Af minimizes 

W2(Af) - ^Af -^-^Af, (16) 
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subject to the constraint that it should be statically ad- 
missible with load increment AF°^*. Minimal values in 
(fT5| and (fT6|) are opposite to each other, and are identi- 
fied in the limit of large systems with the corresponding 
macroscopic elastic energy, i.e. 



^ A 



in which C denotes the 4th rank tensor of elastic moduli. 

Those variational properties are analogous to classical re- 
sults in elasticity of heterogeneous continua [111 , and will 
be used in Section |IV] (technical details being provided 
in Appendix [DJ. 



III. PROPERTIES OF EQUILIBRIUM 
CONFIGURATIONS 

We summarize here some necessary information about 
configurations assembled by different methods, as re- 
ported in paper I [ij, and then isotropically compressed 
to various level of pressure, as described in paper II 0. 
Useful notations and properties are also briefly recalled. 



A. Sample preparation and compression 

Four different configuration series which, as in papers 
I and II, we keep referring to as A to D, were prepared 
under a rather low pressure (k — 39000, corresponding 
to 10 kPa for glass beads, or k = 181000, correspond- 
ing to 1 kPa), and then quasistatically compressed up to 
100 MPa (k ~ 80), with friction coefficient /i — 0.3 in the 
contacts. 

They are characterized in terms of solid fraction $, co- 
ordination number z* of the force-carrying structure or 
backbone of the packing, proportion of rattlers xq (those 
grains do not participate in force transmission at equi- 
librium), normal force distribution, friction mobilization, 
and geometric data such as distribution of interneighbor 
gaps and local order parameters, z* relates to the global 
coordination number z = H^^jn as z* = z/{l — xq). 

Configurations A, B and D were assembled on com- 
pressing a granular gas. Configurations C are obtained 
from A and are supposed to mimic, in a simplified way, 
the dense states obtained by vibration. 

Type A samples are assembled without friction, and 
correspond to the ideal "random close packing" state 
(RCP), which according to the available numerical ev- 
idence is uniquely defined, provided the compaction pro- 
cess is fast enough to avoid all incipient crystalline or- 
der nucleation [1, Section III]. Their solid fraction, ac- 
cordingly, is slightly below 0.64 at low pressure, while 
the coordination number is close to 6, with few rattlers 
{xq ~ 1.5%). Type A configurations may thus be re- 
garded as a simple model for grains that are perfectly 
lubricated in the assembling stage, but such that dry in- 
tergranular contacts have a frictional behavior = 0.3) 



in quasistatic compression. As a variant, another set of 
samples, which we denoted as the AO series, was prepared 
on quasistatically compressing the solid samples without 
friction. 

B states are similar to A ones, except that they are as- 
sembled with a small coefficient of friction, fiQ — 0.02, as 
a crude model for imperfect lubrication in the fabrication 
stage. B states have a smaller solid fraction, $ ~ 0.625 
instead of 0.637 for k < 10^^, and a slightly smaller co- 
ordination number: z* ~ 5.8 instead of 6 at low pressure. 
D states are the loosest of the four series, with $ ~ 0.593 
under 1 kPa, less contacts {z* ~ 4.55 at low P) and more 
than 10% of rattlers. Remarkably, the density of vibrated 
C states is close to the RCP value ($ ~ 0.635), but their 
contact networks are as tenuous as D ones (z* ~ 4.55 at 
low P), with even more rattlers (13%) at low pressure. 
C configurations are thus denser than B ones, but much 
less coordinated. 

The evolution of z* and xq in a pressure cycle up to 
100 MPa, and then back to the initial value, are studied in 
paper II 0, Figs. 1 and 2]. While density increases with 
P, so do the coordination numbers, most notably above 
a few MPa, but upon decompressing many contacts are 
lost and coordination numbers, if initially high, as in the 
A and B cases, end up with much lower values, similar 
to those of poorly coordinated C and D states. Solid 
fraction <I> displays very little hysteresis in such pressure 
cycles. AO (frictionless) states behave very similarly to 
A ones in compression, but do not lose their high coordi- 
nation number on reducing the pressure. The reader can 
refer to papers I and II for more details {e.g., information 
on force distribution and friction mobilization). 



B. Moments of normal force distributions 

In paper I, reduced moments of the distribution of nor- 
mal forces TV were introduced, with the following nota- 
tion: 



Z{a) 



(TV)"' 



(17) 



Another quantity, closely related to Z(5/3) is useful to 
evaluate elastic energies from contact forces. If rxN is 



the ratio 



N 



in any contact, then we define 



Z(5/3) 



5rt 



(iV)5/3 



(18) 



ax denotes the stiffness ratio defined in (|4]). 

It is convenient to use those definitions to write down 
averages of various quantities associated with the con- 
tacts and proportional to some power of the normal force. 
As a useful example, let us relate the average normal stiff- 
ness of contacts to the pressure. This average, from ([3]), 
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reads 
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The sum, running over all contacting pairs is then 
transformed, using Eqns [21 [I7| (with a = 1/3) and ex- 
pressing the contact density as 

Nc _ 3z^ 

into a relation conveniently displaying the dependence on 
pressure, solid fraction and coordination number: 



31/3^ 



7ri/3api/3 



(19) 



This expression of (Kn) will be used on estimating elastic 
moduli. 



C. Degree of force indeterminacy 

In paper I, we also discussed whether the degree of 
force indeterminacy of equilibrated packings could ap- 
proach zero in frictional packings in the rigid, k — > -l-oo 
limit - this being a known property of frictionless sys- 
tems 32]. While the degree of force indeterminacy H is 
directly related to the backbone coordination number z* 
in frictionless packings, for which 



H= in(l-a;o)(z*-6). 



(20) 



its value is more exactly evaluated, for non-vanishing in- 
tergranular friction coefficients, on defining a slightly cor- 
rected value of z*, denoted as z**: 



z* + 



2X2 



(21) 



where X2 is the proportion of 2-coordinated grains (which 
are involved in the mechanism motions mentioned in 
Sec. Ill CI and Appendix [B|. Then the degree of force 
indeterminacy is given by 



A. Numerical results 

Elastic moduli of equilibrated configurations are eval- 
uated as indicated in Section III C 21 Each data point 
on the graphs, throughout the sequel, is based on sev- 
eral macroscopically equivalent load vectors in all 5 avail- 
able samples for each one of the investigated macroscopic 
states. Fig. [T] displays on logarithmic plots the pres- 
sure dependence of shear and bulk moduli in all series 
A, AO, B, C and D during the first compression. Fig. [1] 




1000 



10" 
P (kPa) 
(a)B versus P 
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(b)G versus P 



10* 



10° 



H= -n(l-xo)(z* 



4). 



(22) 



X2 values raise to about 2.5% in configurations C and 
D, in which the ratio of the degree of force indetermi- 
nacy to the number of backbone degrees of freedom, 
i.e., H/[6n(l — xo)], does not decrease below 14%. Only 
in packings assembled with an infinite friction coeffi- 
cient P, [3^ , called Z configurations in paper I, did we 
obtain nearly vanishing H values (H/[6n(l — a;o)] decreas- 
ing to about 3/100). 



FIG. 1: (Color online) Bulk modulus B (a) and shear modulus 
G (b) versus confining pressure P for series A (red, crosses, 
continuous line), AO (red, round dots, dotted line), B (blue, 
asterisks, dotted line), C (black, square dots, continuous line) 
and D (green, open squares, dotted line). Note that results 
for A, AO and B are hardly distinguishable. The dashed blue 
line marked "KJ" corresponds to some experimental data 
between 50 and 400 kPa commented in Section rvD 

clearly shows that the moduli are primarily sensitive to 
coordination number, with well coordinated samples A, 
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B (and AO) displaying larger moduli than C and D, in 
which the contact network is more tenuous {z* ~ 4.5 
and z ~ 4 under low pressure Moduh are much 

less sensitive to packing fraction $: C and D results are 
close to each other at low pressure, when — 0.638 
and $D ^ 0.594 P, They are not strongly influ- 
enced either by the width of the force distribution: A 
and AO states have almost the same moduli (only some 
values of G below 100 kPa differ by more than 5%), 
whereas the probability distribution function of normal 
forces strongly differ as pressure grows (see paper II [I, 
Fig. 5]). 

The increase of elastic constants with pressure natu- 
rally stems from the dependence of contact stiffnesses 
on the force they transmit, as expressed by Eqns. ([3]) 
and ([3]), and due to relation ^ the typical contact stiff- 
ness grows as P^^^ (see Fan. [T9|) . which is the expected 
pressure dependence for macroscopic moduli. Power 
laws are often used to relate elastic moduli to confin- 
ing stresses [l^, I21I, [2^ |4^ , and possible origins for the 
observation of exponents larger than 1 /3 (as on Fig. [T|) 
have been discussed by several authors [40, |41| . One pos- 
sible explanation is the creation of new contacts under 
the effect of the increase of the confining pressure, which 
leads to a denser, stiffer contact network. This mech- 
anism appears in particular to account for the pressure 
dependence of elastic moduli in regular, crystal-like ar- 
rays of identical particles as in the experiments described 
in Refs. and [J. Because of the slight lattice distor- 
sions obtained with imperfect and slightly polydisperse 
spheres, the contact coordination number, which is lim- 
ited, in the rigid limit of k cx3, to 4 in 2D and 6 in 
3D [3^ , is smaller than the nearest neighbor coordination 
number on the dense lattices studied (such as 12 for FCC 
in 3D 5] and 6 for the 2D triangular lattice Q). This 
leaves a large number of neighbor pairs at a distance re- 
lated to the width of the particle size distribution, where 
additional contacts are induced by higher pressures. This 
has been shown by numerical simulations ^42j to produce 
a pressure dependence of moduli closer to P^/^ in some 
pressure range, a phenomenon predicted in part by a the- 
ory presented in [isj . With general, amorphous pack- 
ings, the situation is different because distances between 
neighbors that are not in contact are no longer related 
to a small polydispersity parameter, but are distributed, 
approximately as a power law in some range (see paper 
I [1[), in a way that is characteristic of the disordered 
geometry. Departures from the P^/^ scaling are larger in 
low z states (Fig. [T|), and the largest in C configurations, 
in which contact gains under growing P are faster than 
in D ones. However, apparent power laws with expo- 
nents larger than 1/3 are observed at very low pressures, 
when, from paper II [2, Fig. 2a], the increase of z with P 
is rather slow. Moreover, in the case of C and D systems, 
the exponent of the power law fit for the pressure depen- 
dence of shear modulus G is significantly larger (about 
0.5) than the one for bulk modulus B (about 0.4). These 
features are discussed in paragraph llV Bl below. Changes 



of ratio G/B as P grows are equivalent to changes of the 
Poisson ratio of the granular material, given by 



W-2G 
6B + 2G' 



(23) 



h'* decreases only slightly as P grows for well coordinated 
states A and B, from v* ~ 0.13 at P=10 kPa to v* ~ 0.09 
under 100 MPa. Its larger variations in poorly coor- 
dinated configurations C and D, for which it decreases 
from 0.3 to about 0.1 in the same range, corresponds to 
G increasing with P faster than B. 



B. Simple prediction schemes and relations to 
microstructure 

The simplest approximation scheme to estimate the 
values of elastic moduli, knowing the density and the co- 
ordination number, is based on the assumption of homo- 
geneous strains (or, equivalently, of affine displacements). 
It was introduced, e.g. in [2^, and it is also used by 
Makse et al. in Refs. [1, [l^l (where it is called an effec- 
tive medium theory). It amounts to evaluating the stress 
increments corresponding to strain e using formula ([T]), 
in which the contact force variations are evaluated, via 
Eqn. 15]) , with relative displacements given by 

- Uj = I • {Vj - Ti). 

Using the isotropy of the distribution of contact orien- 
tations, and replacing all normal forces by their average 
value in the computation of contact stiffnesses, this re- 
sults, using relation ([2]), in the following estimates. 



2 I 37r 



2/3 



pl/3 



(24) 



^ 6 + 9aT ^c 
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One thus finds the expected P^/^ dependence, and ob- 
tains moduli proportional to (z$)^/'^. Formulae also 
predict a constant G/B ratio, and thus a constant Pois- 
son ratio: 



= - 0-032 

26 9aT 



(25) 



This latter estimation is considerably smaller than the 
measured values which are given above (shortly after 
Eqn. [23|l . as noted in [8|. This mainly stems from the 
inaccuracy of the estimated value of G [l2|, as we shall 
see. Eqn. [24] suggests to represent ratios 



B 



^2/3pl/3 

G 



(26) 



^2/3pl/3 

as functions of (z*!')^/'^, which is done on Fig. [2] Fig. [2] 
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FIG. 2: (Color online) Reduced moduli, as defined in (|26|l . 
in units of E'^I^P^I^ as functions of (2$)^/'', same colors and 
symbols as on Fig. [1] The estimates given in (I24|l are plotted 
as straight dotted lines. Moduli are plotted for all configura- 
tions in the pressure cycle, showing an approximate collapse 
on a single curve. 



shows that is a siKnificant ly p oorer estimate of G than 
of B, as already noted in [i2[ , for samples of type A or 
B, and even more so in low coordination number configu- 
rations C and D. It also shows that the elastic moduli, as 
a first approximation, can be thought of as determined by 
z and <I>, the former quantity, as it varies more between 
different sphere packings, being the most influential. The 
present study can thus be regarded as a first step towards 
a method to infer coordination numbers nfrom the elastic 
moduli. As stressed in paper I, coordination numbers are 
virtually inaccessible to direct measurements. (It should 
nevertheless be recalled that the present work is limited 
to isotropic configurations, implying isotropy of fabric as 
well as isotropy of stresses) . In this respect it is interest- 
ing to note that the configurations of lower coordination 
number obtained upon decompressing A and B ones af- 



ter they first reach a high pressure level (see paper II 0) 
yield data points on Fig. [5] that stay close to the C and 
D ones corresponding to the same product z$. 

An interesting alternative to the direct use of for- 
mula IT]) is to exploit the variational property expressed 
by (Uni), as explained in AppendixIHl This shows that B'^ 
and G"^ are upper bounds to the true moduli. Account- 
ing for the distribution of forces (see the derivation of 
Ean.[T9)l. those bounds can be slightly improved, yielding 
the analogs of the Voigt upper bound for the macroscopic 
elastic moduli of a mesoscopically disordered continuous 
material: 



B < ^B^Zil/S) 
G < G^°'s* = G"=Z(l/3), 



(27) 



where Z{l/3), as defined in (fT7|) . is always strictly smaller 
than 1. For the bulk modulus, one can also take advan- 
tage of the second variational property, expressed by (jl6p . 
As explained in Appendix IdJ this requires a trial set of 
contact force increments which balance the applied load 
increment. When the stress increment is proportional to 
the preexisting stress, one may take increments of con- 
tact forces that are also proportional to their initial val- 
ues. No such forces balancing a shear stress are available 
for isotropically prestressed configurations. One thus ob- 
tains (see Appendix [D] for details) a lower bound for B 
which is analogous to the Reuss bound for the macro- 
scopic elastic bulk modulus of a mesoscopically disor- 
dered continuous material. This lower bound B^'^^^^ in- 
volves the dimensionless quantity Z(5/3) defined in (fT^ . 
and enables one to bracket the bulk modulus: 



B" 



Z(5/3) 



B 



Rcuss 



< B < B 



Voigt 



B^Zil/^). (28) 



The ratio of the upper bound to the lower one is there- 
fore related to the shape of the distributions of normal 
forces and to the mobilization of friction. Fig. [3] displays 
ratios B/B", B^°''^^/B^, and B^"'''''' / B^ versus (grow- 
ing) pressure P in configurations A and C. These data 
show that in both cases of high (A) and low (C) initial 
coordination number, the bracketing of B given by (|28p 
is quite accurate, the relative difference between upper 
and lower bounds staying below 10% except at the low- 
est pressure for A, and around 15% for C. The Reuss 
estimate is better than the Voigt one in general. It is 
even excellent in the A case for all but the two lowest 
pressure values studied [B appears to be slightly smaller 
than its lower bound at very high pressure because we 
neglected the reduction of intercenter distances between 
contacting grains in the evaluation of the bounds). This 
estimate (see Appendix ID|) becomes exact when the trial 
force increments used are the right ones, and assumes 
therefore that the shape of the normal force distribution 
and the level of friction mobilization do not change with 
pressure. The observation, reported in paper II that 
such quantities are indeed nearly independent on pressure 
in A states (except at the lowest pressure), is therefore 
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FIG. 3: (Color online) Ratio B/B'^ (symbols connected with 

continuous line, error bars), and its Voigt and Reuss bounds FIG. 4; (Color online) A , as defined in ([29}, versus backbone 

(symbols connected by dotted lines) in configurations A (red, coordmation number z* , for isotropic stress increments (red 

crosses) and C (black, square dots), during compression. B asterisks, green for Z states) and pure deviatoric ones (black 

and D samples respectively behave similarly to A and C. crosses, blue for Z states). 



consistent with the success of the Reuss type estimate for 
B in that case. 

Yet, for G no Reuss estimate is available, and the use of 
the Voigt one with the factor Z(l/3) hardly reduces the 
discrepancy between and G: this factor can be read 
on Fig. [21 where it coincides with the upper bound, and 
hence G'^°'st jg only smaller than G'^ by a few percent. It 
thus overestimates the true shear modulus by 30 to 40% 
in well-coordinated states, and even by a factor of 3 in 
poorly coordinated ones at low pressure. 

C. Fluctuations and more sophisticated prediction 
schemes. 

The Voigt or mean field approach ignores fluctuations 
in grain displacements and rotations. As an indicator 
of the amplitude of such fluctuations is the average of 
squared particle displacements, we measured the ratio 

in which the squared strain amplitude, ||e|p = Ci+^i+^ij 
normalizes the average of squared displacement fluctu- 
ations, as defined in ([7]), the sum running over the 
n* = nil — xq) force-carrying grains. Fig. 3] displays 
the values of evaluated in all samples at the different 
values of the confining pressure. is distributed over 
some fairly wide interval in similar configurations, but 
is systematically larger for purely deviatoric stress incre- 
ments than for isotropic pressure steps and has a clearcut 
decreasing trend as a function of backbone coordination 
number z* . (To add data points with lower z* on Fig. [3] 
we also used the Z series, infinite friction samples, as 



described in paper I and Section IIIIl isotropically com- 
pressed and equilibrated at five pressure levels, from 1 to 
100 kPa). This is consistent with the approximation that 
ignores fluctuations being less accurate for shear stresses 
than for isotropic ones. 

More elaborate prediction methods for elastic moduli 
were proposed. Kruyt and Rothenburg [53, SI] consid- 
ered two-dimensional assemblies of non-rotating parti- 
cles, and succeeded in applying a variational approach 
such as the one we use for bulk modulus B to the eval- 
uation of shear moduli as well. Velicky and Caroli [i^ 
studied the case of an imperfect lattice system with con- 
tact disorder, as in the experiments of ^5i| and the simu- 
lations of [13] ■ Jenkins et al. [2^ dealt with frictionless 
sphere packings. More recently. La Ragione and Jenk- 
ins [45] published an approximation scheme which is di- 
rectly comparable to our simulation results, the results 
of which are denoted as LRJ below. 

Those estimation procedures improve upon the Voigt 
assumption that relative displacements are ruled by the 
average strain on considering small sets of displacements 
and rotations, either associated to one grain, or to a con- 
tacting pair. Those degrees of freedom are allowed to 
fluctuate while their surroundings abide by the Voigt as- 
sumption. Optimal values of the fluctuating variables 
are then to be determined on solving the corresponding 
system of equilibrium equations for the selected small 
set of degrees of freedom. Such approaches necessarily 
involve complex treatments of the random geometry of 
local grain arrangements, especially on attempting to ex- 
press the predicted moduli with a limited amount of sta- 
tistical data. In [4^ we numerically check some of the 
approximations involved, on exactly solving the required 
set of local equilibrium problems. We show 46] that the 
discrepancy between observed and predicted shear mod- 
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uli is reduced, down to 50% in the worst cases of C and 
D samples under low pressures. 



coordinated systems, for which we now test another kind 
of theoretical prediction. 
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FIG. 5: (Color online) Ratio of estimated to measured values 
of shear moduli, G"="VG, with G°=' = G^°'«^' (larger values, 
dots connected by dashed lines) and C^^* — G^^^ (smaller 
values, dots connected by solid lines), versus P or , in 
sample series A (red), B (blue), C (black) and D (green). 

Leaving aside the discussion of the various approxima- 
tions involved in the derivation of the LRJ predictions 
(some important steps of which are tested in [i^), we 
now confront them with the numerical data. 

We observe that the LRJ formulae do not improve 
the predictions of bulk moduli over the Voigt and Reuss 
bounds ((28|) . Occasionally, they even predict values below 
the lower bound. Yet, as shown on Fig. [5l the estimates 
QhRj obtained for shear moduli are much better than the 
Voigt ones (|27p . Shear moduli in well-coordinated states 
A and B are accurately predicted, while the discrepancy 
in poorly coordinated systems C and D, from a factor of 
three with the Voigt formula, are down to about 50%- 
70% with the LRJ one under the lowest pressure levels. 

LRJ formulae yield moduli that are proportional to av- 
erage contact stiffnesses (in which we introduced an ad- 
ditional factor of Z{l/3), see Eqn. [T^ . with coefficients 
involving rational functions of the backbone coordination 
number z* , and also the variance of the fluctuations in 
the number of contacts of backbone grains. The LRJ pre- 
dictions still overestimates the shear modulus in poorly 



V. THE CASE OF POORLY COORDINATED 
NETWORKS 

The specific elastic properties of configuration series C 
and D, with their small coordination numbers, are rem- 
iniscent of frictionless packings ^, '47*1, in which a sim- 
ilar anomalous behavior of G as a function of pressure 
has been reported. Here we review these properties of 
packings with no tangential forces (Section I VA[) . and we 
discuss a possible explanation (48j . and its applicability 
to poorly coordinated frictional packings (Section IV B|) . 



A. Frictionless packings 

Although samples of series AO were confined with no 
mobilization of friction, elastic moduli shown on Fig. [1] 
have been computed with taiigential elasticity in the con- 
tacts, just like, e.g. in Ref. [l3|. It is assumed for state 
AO that friction is not mobilized in the preparation pro- 
cess. In other words samples are perfectly annealed to 
a local minimum of mechanical energy inconfiguration 
space, but the response to some stress increment implies 
tangential forces in the contacts. Results are of course 
different if contacts are still regarded as frictionless on 
evaluating elastic properties. Fig. [6] compares this new 
set of values, which we denote as AOO, to AO ones. Bulk 



moduli (Fig. 6(a)) are only slightly higher (about 10% 
at low pressure) with tangential elasticity. AOO values 
correspond to frictionless packings, in which a relation 
between pressure P and the increase in solid fraction A<i> 
above the rigidity threshold {i.e., above the value of <i> 
in the limit of vanishing pressure) was obtained in paper 
I [T, Eqn. 31], as a direct consequence of the isostaticity 
property of the backbone. It is straightforward to check 
that this relation, on taking the derivative with respect 
to <& and relating modulus B to P, yields the Reuss type 
estimate of -B, as expressed by ([25)1 and in which z 
and $ are replaced by their values in the limit of P — > 
(which does not significantly affect the result at low pres- 
sure) . Both approaches are based on the assumption that 
forces vary proportionnaly to P, which is a particular 
consequence of isostaticity, hence the exactness of the 
Reuss estimate for AOO results. The distribution of force 
values in sample series AO, once normalized by the ap- 
plied pressure level, was checked in paper II to remain 
very nearly constant. 

The small influence of tangential elasticity on bulk 
moduli, which is responsible for the difference in B val- 
ues between AO and AOO series, is not surprising, as both 
the Voigt and the Reuss- like approaches, which restrict 
the values of B to the interval given by lead to the 
assumption of vanishing tangential forces. 

The shear modulus, on the other hand, as noted in 
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FIG. 6: (Color online) Elastic moduli of samples AO prepared 
at different pressures without friction, computed witli (AO) 
and without (AOO) tangential elasticity, AOO results corre- 
sponding to completely frictionless packings. 



refs. [8|, |47|, is singular in frictionless packings under 
isotropic stresses: while values of G in the AO series 
vary approximately proportionally to B, and are of the 
same order, as observed above in Sec. lIV A] shear moduli 
of frictionless systems AOO (Fig. |6(b) ) are considerably 
smaller, and vary faster with P. This increase is very 
well fitted by a power law with exponent 2/3, in agree- 
ment with [43] ■ An explanation for the singular behavior 
of G is suggested in [1^, |4^, as follows. First, some of 
the pressure dependence of G is simply due to the influ- 
ence of pressure on average contact stiffness (-ft' at). One 
should therefore rather explain the pressure dependence 
of G/{Kn)- Then it is argued that this amplitude is 
proportional to the degree of force indeterminacy, or to 
z* — 6. More precisely, the shear modulus should scale 
as the degree of force indeterminacy per unit volume, or 
equivalently as {z* — 6)(1 — x^)^. This is the crucial 
part of the argument. Leaving aside a discussion of its 




FIG. 7: Reduced shear modulus Qa, defined by l|30p . in 
frictionless configurations AOO versus backbone coordination 
number 2*. The straight line is the best linear fit through the 
6 leftmost data points. 



justification (which would require detailed calculations 
of sets of self-balanced contact forces and response func- 
tions within the contact networks) we check here for its 
practical validity. To do so, we define a reduced shear 
modulus ga on dividing by (1 — a;o)$ and by (Kn)- We 
thus have, from (fT9l). 



5a 



i^2/3pl/3^(l/3)(l_^g)$2/3' 



(30) 



and we test (Fig. [7]) whether it varies linearly with z*. 

The final part of the demonstration of Ref. (see 
also nil), suggests to evaluate the increase of coordina- 
tion number with pressure on relating both quantities to 
the increase in packing fraction above rigidity threshold 
A$. Such a relation between P and A<i> in isostatic fric- 
tionless packings was written in 1], with P cx (A<1>)'^/^. 
The additional ingredient is a scaling form of the incre- 
ment of z* with A<i>: 



6(x A$i/^ 



(31) 



A homogeneous shrinking approximation is suggested 
in [4^ to derive relation (|3ip . based on the assumption 
that the gap-dependent near neighbor coordination num- 
ber z{h) grows like z{h) — 6 oc h'^-^. However, as shown 
in [l| we observed an exponent 0.6 instead, and our data 
therefore do not confirm this argument. 

Nevertheless, the proportionality of the singular ampli- 
tude ga of the shear modulus to z* — 6 is accurately satis- 
fied, as shown on Fig.[71 The linear fit of the dependence 
of ga on z*, through the 6 first data points, is very good 
and predicts a vanishing modulus for z* = 5.994 ± 0.008. 
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B. Packings with intergranular friction 

In paper I we concluded that frictional packings pre- 
pared in low coordination states did not approach iso- 
staticity under low pressure. However, one may test 
whether the amplitude Qa varies linearly with the degree 
of force indeterminacy when it is small enough, even if 
it does not approach zero. The Z states, on the other 
hand (see Section |TTT] and paper I [l|), were prepared 
with an infinite friction coefficient and have nearly van- 
ishing force indeterminacy at low pressure. Fig. [51 in 
which Qa data for states C, D and Z are plotted versus 
the corrected backbone coordination number z**, which 
determines the degree of hyperstaticity per degree of free- 
dom by (I22p . shows that the linearity is very well sat- 
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FIG. 8: (Color online) Amplitude ga (Egn. l30|l versus z** , as 
defined in (|22p . in samples of types C (black), D (green) and 
Z (5 lower data points, red). The dotted line is a linear fit to 
the 4 lowest Z data points. The LRJ predictions are shown 
for series z, C and D, as indicated (asterisks joined by dotted 
lines, same color code). 

isfied for Z states. Z, C and D points lie approximately 
on the same curve, showing that the macroscopic modu- 
lus is controlled by z**, in spite of the differences in the 
structures of states C, D and Z. (Due to the greater mi- 
cromechanical changes observed upon reducing the pres- 
sure from its maximum value , the corresponding data 
points for C and D states lie on a different curve, not 
shown on the figure.) The linear fit still approximately 
applies to the lowest values for D and C. As expected, 
the linear fit predicts G = for z** = 3.99 ±0.02, i.e., a 
shear modulus vanishing proportionnally to the degree of 
force indeterminacy. Z configurations have a larger pop- 
ulation of rattlers and divalent grains than C or D ones: 
xq ~ 0.184 and X2 — 0.068 at 1 kPa. Consequently, 
on fitting ga versus z* instead of 2**, ga would appear 
to vanish unambiguously below z* — 4, for the value 
z* — 3.93 ± 0.02. The correction due to 2-coordinated 
beads improves the accuracy of the result that G van- 
ishes with H . 



Fig. [8] also displays the values of amplitudes ga pre- 
dicted by the LRJ formulae [i^, as discussed in Sec- 
tion [IVC] As should normally be expected for such an 
estimation procedure, based on the local equilibrium of 
one pair of grains embedded in an elastic medium, the 
LRJ approach is unable to capture the vanishing of shear 
moduli in the limit of z** — > 4 (i.e., H 0), since the 
rigidity properties of tenuous networks are determined by 
more collective effects. The low level of force indetermi- 
nacy provides a complementary approach to the estima- 
tion schemes evoked in Scction llV CI to predict the values 
of shear moduli in isotropically compressed packings. 

We conclude therefore that the proximity of a state 
devoid of force indeterminacy, however unreachable, ex- 
plains the anomalously fast increase of the shear modulus 
with the pressure for low coordination frictional pack- 
ings, as observed on Figs. [T] and [2] As to the increase 
of the degree of hyperstaticity, or of z** , with pressure, 
its prediction seems to be even more difficult than in the 
frictionless case. What would be needed is an accurate 
prediction of small changes in z* , which, as observed in Q 
(paper II) the simple homogeneous shrinking assumption 
does not provide, due, in particular, but not only, to its 
inability to deal with the recruitment of rattlers by the 
growing backbone. 

The proximity of a "critical" value of the number of 
contacts on the backbone also entails specific properties 
of the eigenvalues of the stiffness matrix (the "density 
of states" in the language of solid-state physics), with 
a large excess of soft modes [U [13]. A similar behav- 
ior, both for the cigenmodes of the stiffness matrix and 
for some shear elastic constant was observed in [s^ in 
2D simulations of anisotropic states. From this partic- 
ular set of results in anisotropic packings and from the 
Reuss approach to estimate the bulk modulus one may 
deduce that the non-singularity of i?, as opposed to G, 
directly stems from the isotropic state of stress on which 
load increments are applied. On increasing P, in a good 
approximation (the better the lower the degree of force 
indeterminacy), one just rescales the contact force values. 
Load increments that are not proportional to the preex- 
isting load, on the other hand, tend to produce large 
fiuctuations (see Fig. |4]) and soft responses in poorly co- 
ordinated contact networks. 



VI. COMPARISON WITH EXPERIMENTAL 
RESULTS 

The assembling procedures of states B and C (Sec. lIIIj) 
can be regarded as idealized models for lubrication and 
vibration. Jia and Mills [l^, [13] measured sound wave 
velocities in glass bead packings, some samples being 
densified by repeated taps on the container, and oth- 
ers mixed with a very small quantity of a lubricant (tri- 
oleine) . The beads were placed in a cylindrical container 
and then compressed by a piston, transmitting a con- 
fining pressure. Velocities of longitudinal and transverse 
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sound waves propagating in the vertical direction (or- 
thogonal to the piston) were measured in the 70 kPa- 
800 kPa range. Those velocities, which we denote as 
usual as Vp and Vs, relate for an isotropic material to 
bulk and shear moduli and mass density p as 



Vp 



IB 



4 

and Vs 




Fig. [9] displays the sound velocities in both types of sam- 
ples, along with the results of Domenico [si], measured on 
a glass bead sample in a higher pressure range. Within 
the range of vertical pressure P investigated in the ex- 
periments of [13] , the packing fraction of dry beads varies 
between 0.633 and 0.637, while lubricated ones are less 
dense, $ ranging between 0.613 and 0.617. Sound veloc- 
ities are nevertheless larger in lubricated packings. 

Comparisons with our numerical data, also shown on 
Fig. [9l in spite of the differences in preparation proce- 
dures (which are idealized, and involve somewhat arbi- 
trary choices of parameters in simulations) and loading 
(oedometric compression in experiments, isotropic com- 
pression in simulations), reveal some interesting qual- 
itative convergences and semi-quantitative agreements. 
Specifically, we note that: 

• Numerical "lubricated" samples B are also less 
dense, but stiffer than numerical "vibrated" sam- 
ples C. 

• Sound velocities in B samples increase with P 
slower than in C ones, and sound velocities in lu- 
bricated laboratory samples increase slower than in 
vibrated ones. 

• Numerical C samples are better models for dry ex- 
perimental packings assembled by vibration than 
A ones: values of elastic moduli are in much better 
agreement. 

One may therefore attribute the difference in sound ve- 
locities reported in [22] between dry and lubricated pack- 
ings to a difference in coordination number, like in nu- 
merical states B and C. (Su ch an interpretation differs 
from the one set forth in [2^, which relies on the filling 
of open interstices by the lubricant) . 

The traditional numerical route to obtain dense sam- 
ples, i.e. the use of a vanishing or low friction coef- 
ficient as for systems A and B, fails to reproduce the 
elastic properties of dense samples assembled by vibra- 
tion. Those appear to be better simulated with the newly 
introduced numerical procedure resulting in C samples, 
which have a much lower coordination number for the 
same density. Laboratory samples with a solid frac- 
tion approaching the RCP value might well, especially 
if their preparation involves vibrations or tapping, pos- 
sess as small a density of force-carrying contacts as our 
numerical samples of type C (^ ~ 4.05). 



Such a conclusion, in favor of low-coordinated numeri- 
cal samples as better models for experimental dense pack- 
ings of dry beads than conventional, A-type ones, appears 
to contradict the results of Makse et al. [I,[l3]- Those au- 
thors simulated what we denoted as the AO sample series, 
and reported good agreements with Domenico's experi- 
mental results [oj and with their own measurements. We 
checked that the agreement between their numerical re- 
sults and our AO data was excellent. We attribute the 
conflicting conclusions to their comparison being done in 
a much higher pressure range than the one of Jia and 
Mills' experiments: as apparent on Fig. [SI the confining 
pressures in Domenico's experiments are all above 2 MPa. 
Likewise, P values all exceed several MPa in the exper- 
iments performed by Makse et at. In this range, differ- 
ences between A and C samples, as apparent on Fig.[5]for 
sound velocities, as well as on [2j, Fig. 2a] for coordination 
numbers, tend to dwindle as P increases. The discrep- 
ancy between numerical results on A-type systems and 
experimental results on dry bead packings is much lower 
under high pressure. Yet, numerical samples of type C 
still fit the experimental data better. The apparent ex- 
ponent in a power-law increase of sound velocities with 
P, i.e. the slope on Fig. [9l is, in particular, better re- 
produced by C data than by A ones. On discussing such 
a high pressure range, one should nevertheless keep in 
mind the possible occurrence of non-elastic behavior in 
the contacts, as pointed out in (paper II), where the 
maximum stress levels in contact regions were estimated. 

The fast increase of G as a function of P in C samples 
is not observed in the experimental results of Jia and 
Mills. Furthermore, using formula ([^5]) (which assumes 
isotropy), these results correspond to Poisson ratios (be- 
tween 0.32 and 0.34), which are larger than numerical 
results and do not vary with P. However, these data are 
bound to be affected by stress anisotropy. 

In this respect, a comparison of numerical results with 
the data of Kuwano and Jardine @ is easier, as those 
were measured in glass bead samples under isotropic 
stress states, from about 50 to 400 kPa. The sam- 
ples of Kuwano and Jardine have similar densities to 
D ones ($ ~ 0.59), and are initially made by "air plu- 
viation" , i.e., deposited under gravity, in a controlled 
procedure ensuring homogeneity. The values of shear 
moduli are close to the numerical values for C and D 
states, and vary with P even faster, as G oc P'^-^^. The 
power law fit through these data correspond to the line 
marked "KJ" on Fig. 1(b) Kuwano and Jardine, combin- 



ing static small-strain tests and sound velocity measure- 
ments, could evaluate the 5 independent elastic moduli 
of the transverse isotropic granular material assembled 
under gravity. To compare our numerical results, ob- 
tained in isotropically assembled systems, with theirs, 
we ignored the moderate effect of the fabric anisotropy 
on experimental elastic moduli and used, in Fig. l(b)[ 
the moduli corresponding to a shear strain in the ver- 
tical plane, the shear modulus in the horizontal plane 
being about 7% larger. Another similarity between our 
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FIG. 9: (Color online) Sound velocities Vp (upper date points) and Vs (lower data points) versus confining pressure P on 
double logarithmic plot for the experimental dry (vibrated) and lubricated samples of Jia and Mills, and of Domenico, to which 
numerical values for simulated states A, B, C and D are compared. 



results in D or C states and the data of [6] is the pres- 
sure dependence of the two Poisson ratios v^h and Vhh 
which couple stress and strain components in 2 differ- 
ent directions respectively defining vertical and horizon- 
tal planes. Despite some scatter in the measured values 
of these ratios, v^h and Uhh show a marked decreasing 
trend between P=80 kPa and P=400 kPa. Finally, we 
compare the Young moduli mesured in [6*] in vertical and 
horizontal directions to our numerical values in states A 
to D on Fig. [TOl Our numerical values for Young mod- 
ulus E* in systems with low coordination numbers are 
similar to those results, but the pressure variation seems 
faster (with exponent ~ 0.6) in experiments. 



We conclude therefore that, although more systematic 
confrontations with experimental results are necessary, 
some features of the moduli in low coordination numeri- 
cal packings are apparently observed in the rather loose 
glass bead samples of Kuwano and Jardine. 
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FIG. 10: (Color online) Young modulus E* in numerical sam- 
ples A, AO, B, C, D, labelled with same colors and symbols as 
on Fig. |9l as a function of P, compared to fits through data 
points of Kuwano and Jardine (dotted lines, KJ) with two 
sets of values because of fabric anisotropy. 
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VII. ELASTIC REGIME 

The elastic moduli express the relationship between 
small stress and strain increments. We now wish to eval- 
uate the elastic range, and to explore the origins of the 
breakdown of elasticity for larger increments. 

A. Limit of linear regime 

Motivated by comparisons with experiments, we tested 
the predictions of linear elasticity, as evaluated with the 
moduli obtained from the stiffness matrix, versus the 
full MD simulation for small and slowly applied load 
increments, in the case of a triaxial axisymmetric com- 
pression. This test, familiar in geomechanical engineer- 
ing [sil . [s^ . [ssj , consists in increasing one stress compo- 
nent, say CTi, while the other 2 are kept constant at the 
initial value of the isotropic pressure in the initial state: 

= f 3 = P- More often, in the laboratory, one controls 
strain component ei (called "axial strain" ) with the mo- 
tion of a piston, imposing a constant, slow strain rate ei, 
while the lateral stresses are maintained by the pressure 
of a fluid surrounding the sample, which is wrapped in 
an impervious membrane. It is customary to express the 
results of such a test with two curves, representing, as 
functions of axial strain, the deviator stress q = ui — 
and the volumetric strain e„ = ei -|- e2 -I- £3, the initial 
isotropic state being chosen as the origin of strains. 

One should have in the quasistatic regime (small 
enough ei), within the linear elastic range, for small 
enough e, 

q = E*^i , ^ 

= (1 - 1v )ei, 

where E* — — and v* (see Eqn. [^5)1 are respec- 

3i? 4- G 

tively the Young modulus and the Poisson ratio of the 
material in the initial state. 

We simulated triaxial compressions for small axial 
strains and compared the resulting deviator stress and 
volumetric strain curves to p2p . To minimize dynamical 
effects in simulations of quasistatic behavior, the iner- 
tia parameter / was kept below lO^"* or even 10^^ for 
the most fragile, low pressure samples. (/ is defined by 
/ m 

/ = ew — with m denoting the mass of one grain, see 
V aP 

papers I and II P, [I])- Fig. [Til shows the typical results 
of such a comparison in the case of one C sample, un- 
der isotropic initial pressures P growing from 10 kPa to 
1 MPa. This comparison first shows that full MD com- 
putation results do admit, in excellent appproximation, 
as slopes of the tangents to the curves at the origin, the 
appropriate values deduced from the evalution of elas- 
tic moduli by the stiffness matrix approach, as expressed 
by ([5^ . This confirms the statements made in Sec. Ill CI 
and checked in Appendix Elabout the definition of elastic 




(a)Deviator stress normalized by P = 0-3 versus axial strain ei, for 
the five P values indicated. 




(b) Volumetric strain — eu, showing contractance, versus axial 
strain, same P values, growing according to arrow. 

FIG. 11: Deviator stress (a) and volumetric strain (b) versus 
axial strain in beginning of triaxial compression of a type C 
sample. Dots show MD results of triaxial compression while 
straight lines have slopes given by (I32p . The volumetric strain 
is shown with an axis oriented downwards to better visualize 
the decrease in sample volume, which contrasts with its dila- 
tancy for larger strains. 

moduli: this definition makes sense as a very good ap- 
proximation in spite of the slight directional dependence 
of contact stiffness matrix 

Values of t\ for which the q(t\) curve significantly de- 
viates from its tangent at the origin increases from about 
- 10~5 to ~ 10-"* as P grows from 10 kPa to 1 MPa. This 
linear elastic range limit, if expressed with stress ratio 
(7/(73, shows but a slow increase with uz = P, starting 
around (7/(73 ~ 0.05. The initial slope of the volumet- 
ric strain curve increases with pressure, in C samples, in 
agreement with the results on Poisson ratios fSec. IIV X)) . 
Deviations from the linear range of ([5^ appear a little 
sooner, for t\ slightly below 10~^ at 10 kPa, increasing 
to a few times 10"^ at 1 MPa. 
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In the literature on sand properties [1, [Tl| , it is of- 
ten observed that the approximately linear elastic range 
about a prestressed reference state extends to strains of 
order 10~^ or 10~^. On Fig. [TH we plotted the value 
of ei for which the deviator stress starts to deviate from 
P2|) by more than 5%, versus the confining pressure, for 
series A, C, and D. Recalling that most experimental re- 




P (kPa) 



FIG. 12: (Color online) Threshold e?''" above which q differs 
from E*e\ by more than 5%, versus P, for series A (red), C 
(black) and D (green). The dashed line has slope 2/3. Data 
points corresponding to one sample series are connected by 
thin dotted lines. 

suits are collected in the range 50 kPa< P < 1 MPa, 
these data confirm the experimental observations (isj 
on sands in terms of order of magnitudes for all three 
sample series. However, they also witness a systematic 
growth of the elastic threshold ef'^^ with P, roughly as 
P"^/^. This suggests a constant elastic deviator interval 
relative to the confining pressure, q^^'^^/P, on assuming 
E oc P^/^. Figs. [Ql and [HI show that the elastic range 
is better expressed in that form, as the threshold ratio 
gcias^p displays much smaller variations as a function of 
P: unlike Fig.[T21 those graphs show stress intervals on a 
linear scale. Expressing the extension of the linear elas- 
tic regime in terms of strains, as done in the literature 
on sands, allows one to gather the different sample series 
within the same range of magnitudes around 10~^, pro- 
vided the confining pressure stays within the interval that 
is most often investigated in experiments. However, the 
pressure dependence is better accounted for on express- 
ing the upper limit of the linear elastic range in terms of 
stress increments, relative to the confining stresses. The 
trend in low-coordinated systems C and D is an increase 
of the linear elastic interval, expressed as a stress ratio, 
with P. At the lowest pressure, the smallness of this in- 
terval (typically about lO^'^ for ratio q/P in D samples 
under P = 1 kPa) is characteristic of the greater fragility 
of tenuous networks, a phenomenon, once again, remi- 
niscent of the situation of nearly isostatic force-carrying 
structures in frictionless packings under low pressure. In 
the limit of rigid grains, frictionless systems (if their prop- 
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FIG. 13: Stress ratio 5°'^7P above which q/E* differs from 
ei by more than 5%, versus P, for series A. 
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FIG. 14: (Color online) Stress ratio q'^^'^yP above which q/E* 
differs from ei by more than 5%, versus P, for series C (black, 
lower data points) and D (green) . 



erties do not qualitatively depend on space dimension) 
should behave as described in [53 |: in packings of rigid 
disks the deviator stress increments causing exactly iso- 
static 2D contact networks to fail and rearrange were 
found to approach zero as an inverse power of the num- 
ber of particles in the limit of large systems. 

Series A configurations, on the other hand, have de- 
creasing linear elastic intervals as a function of pressure. 
This is due to the increase of friction mobilization, which 
leads to larger non-elastic terms in the response to load 
increments. 



B. Onset of irreversibility 

So far we have been testing the accuracy of linear elas- 
ticity, i.e., the predictions of p2p with the initial mod- 
uli. It is known from experiments that granular mate- 
rials cease to be elastic outside this linear regime. This 
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FIG. 15: (Color online) Deviator stress (black, left axis) and 
volumetric strain (red, right axis) versus axial strain in triax- 
ial compression for small strains, with unloading curves, in a 
D sample initially under P=100 kPa. The blue dotted curves 
show the results of a calculation with the sole contacts that 
are initially present. 



may be checked on testing for reversibility in a strain 
cycle. The effects of unloading from various points on 
the triaxial compression curve are shown on Fig. 1151 in 
a type D sample with (T2 = = 100 kPa. On revers- 
ing the sign of ei, while still maintaining constant lat- 
eral stresses a2 — crs, one observes that the unloading 
curve starts with a slope close to the initial slope of the 
loading curve. This common slope is the elastic mod- 
uli corresponding to the linear response for very small 
strain increments. Consequently, the response to a devi- 
ator stress is no longer reversible as soon as it ceases to 
be linear, as shown on Fig. I15i on which the axial strain 
that results from a deviator stress cycle is represented. A 
large proportion of the strain entailed by a relative stress 
increase of order 10"-'^ is not recovered on returning to the 
initial load. The response to an isotropic load increment 
is much closer to reversibility, even for much larger stress 
variations, as shown on Fig. 1161 which displays the stress- 
strain curve in an isotropic pressure cycle. As P varies by 
a factor of 2, about 93% of the volume increase is recov- 
erable upon unloading to the initial pressure. Only for 
the large pressure cycles as studied in paper II ^2., Fig. 2] 
can one observe notable irreversible changes, in coordina- 
tion number rather than in density. One thus finds again 
that the response to increments in load intensity (here: 
isotropic compression) strongly differs from the response 
to changes in load direction (here: deviator stress). 




FIG. 16: Isotropic pressure increase versus relative volume 
change in pressure cycle, P growing from 10 to 20 kPa in C 
samples, and then decreasing back to 10 kPa. 



C. Origins of nonlinearity and irreversibility 

For deviatoric stress increments, the gradual onset of 
irreversibility, as the applied load variation increases, co- 
incides with the breakdown of linearity. The lack of re- 
versibility, as shown on Fig. I15[ has two different origins. 
One is the mobilization of friction, and the second is the 
failure of the contact network: the packing eventually 
breaks apart and rearranges. In order to detect the oc- 
currence of this latter kind of event, one may compute 
the response in the beginning of a triaxial compression 
of some samples with MD calculations in which only the 
initially present contacts, in the isotropic state, are taken 
into account. One thus tests the ability of the initial con- 
tact network to support different stress values. One then 
observes, as shown on Fig.[T5l that the initial contact net- 
work proves unable to support a deviator stress beyond 
a certain limit: the q versus ei curve reaches a maximum 
if ei is controlled. This witnesses the propensity of the 
packing to become unstable and gain kinetic energy be- 
fore it finds a new contact network that is able to support 
a larger deviator. This happens the sooner the larger the 
stiffness parameter k in systems with low coordination 
numbers. Thus the corresponding ratio q/cr^ decreases 
from about 0.2 in the example of Fig. [T5l corresponding 
to one D-type sample under 100 kPa, to, typically, 0.1 un- 
der 3 kPa, and 0.05 under 1 kPa. In this respect (as for 
the values and pressure dependence of shear moduli) low 
coordination frictional packings exhibit, in a weakened 
form, similar singular behaviors to frictionless ones. The 
opposite behavior is that of A-type packings, with a very 
large coordination number. As shown on Fig. 1171 the ir- 
reversibilities are smaller and the initial contact network 
proves able to withstand considerably larger relative de- 
viator stresses, q/a^ ^ 1.1 under 100 kPa. The behavior 
of A packings under shear is therefore somewhat inter- 
mediate between that of low coordination systems C or 
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FIG. 17: (Color online) Analog of Fig. 1151 in a sample of type 
A. Note the different scales. 



D under not too low confining stress, and the response 
to isotropic pressure increases. The behavior of A sam- 
ples under triaxial compression is to a large extent de- 
termined by the response of the initial contact network, 
and the rise of deviator g as a function of axial strain is 
very fast. This is a typical feature of well-coordinated 
packings, as studied in [5^. In Refs. [EHl and [s^, two 
different types of strain are distinguished: those due to 
the deformability of contacts, and those stemming from 
network failures and rearrangements. As long as the first 
type of strains is the only one present, the behavior is 
close to that of the inital contact network. Beyond the 
stability range of the initial contact network, the effect of 
rearrangements dominate [s^, and strains are produced 
by local instabilities which can be described with a rigid 
grain model [STj . 

Refs. [Sl, [sll show that well-coordinated isotropic 
packings (2D analogs of A systems) can support rather 
large deviator stress increments in the k —^ oo limit, 
whatever the sample size. Systems with lower coordi- 
nation numbers appear to exhibit intermediate behav- 
iors between this one and the "fragility" , defined as the 
propensity to rearrange for arbitrary small stress incre- 
ments in the large system limit [H, [13, HE HI], of as- 
semblies of rigid, frictionless grains. The stability range 
of given equilibrium contact networks extends to smaller 
stress increment intervals in C or D-like packings, but 
we expect it to remain finite for arbitrary low pressure 
levels. These properties, and the distinction of two types 
of strains, are further explored and discussed in a forth- 
coming paper [60j . 



VIII. CONCLUSIONS AND PERSPECTIVES 

Our numerical results can be summarized as follows. 

Elastic moduli of granular packings are primarily sen- 
sitive to the stress level, via the average contact stiffness, 
which is proportional to P^/'^(z<l>)~^/^ under pressure P. 
Once this effect is taken into account, important differ- 
ences remain between the elastic properties of different 
packing structures, and systems assembled with the same 
density might exhibit large variations in their moduli, 
since those are essentially related to coordination num- 
bers. Under isotropic pressures, one should distinguish 
between bulk and shear moduli. The bulk modulus, in 
all studied configurations, is efficiently evaluated by the 
Voigt and Reuss-like bounds, with an error smaller than 
20%. In general, we expect a difference between the re- 
sponses to changes in load intensity on the one hand, 
and to changes in stress direction on the other hand. In 
isotropic systems, the latter correspond to purely devi- 
atoric stresses, the effect of which differs the most from 
that of hydrostatic pressure increments in low coordina- 
tion contact networks. In such cases shear moduli are 
anomalously small and increase faster with the confining 
pressure. In well coordinated states, such as A and B, 
satisfactory estimates of B and G moduli on using im- 
provements of the Voigt approximation, based on locally 
independent fluctuations about average strains, such as 
the La Ragione- Jenkins (LRJ) approach. Moreover, the 
additional stiffening effect of the increase in coordina- 
tion number, due to compression, might, for pressures 
in the MPa range, be reasonably predicted with the ho- 
mogeneous shrinking assumption, or similar refinements 
thereof. Such schemes nevertheless require rather de- 
tailed statistical knowledge of local particle configura- 
tions. On the other hand, the shear response of low co- 
ordination packings, such as C and D, is better described 
with reference to a state with no force indeterminacy, 
even though such a state is not closely approached in the 
limit of low pressures (or large stiffness level, k —^ +oo) 
except for unphysically large friction coefhcients (as in 
the Z configuration series). In the rigid, k +oo limit, 
shear moduli become proportional to the level of force 
indeterminacy, which directly relates to z* — 4, with a 
small correction due to divalent particles (hence the def- 
inition and use of z**). The dependence of z* — 4 on 
pressure seems however difficult to predict with the nec- 
essary accuracy [i]. The physical origin of the breakdown 
of linear elasticity beyond an interval of very small strains 
depends on z* as well. Linear elasticity fails considerably 
later, in terms of relative stress increments, for isotropic 
pressure changes, and, to a lesser extent, in high coordi- 
nation packings subjected to deviatoric stress paths. In 
such cases linear elasticity breaks down because of fric- 
tional forces. Elasticity ceases to apply for very small 
shear stress increments in low coordination systems, the 
smaller the closer they are to packings devoid of force 
indeterminacy. In such cases the main physical origin 
of nonelastic response is network fragility, as the contact 



19 



structure breaks apart in response to stress increments. 
Extension of linear elastic regimes observed in numeri- 
cal simulations agree semi-quantitatively with observa- 
tions on sands. The shape of the stress-strain curves 
beyond the elastic range correlates with the coordination 
number, with a much stiffer response in well coordinated 
packings. On comparing numerical and experimental re- 
sults, the low pressure regime of poorly coordinated net- 
works corresponds to the lowest pressures for which lab- 
oratory results are available. The characteristic features 
of this regime, such as G increasing with P faster than B, 
or Vs faster than Vp, as apparent in C and D numerical 
configuration series, are not observed in the experimen- 
tal data of Jia and Mills [H, 113] on dense samples. Yet 
some other measurements made by Kuwano and Jardine 
on looser sphere packings Q show similar trends to C 
and D-state simulations. 

The variety of observations corresponding to the same 
pressure and density values for the same material con- 
firms the sensitivity of elastic moduli to otherwise un- 
detectable differences in inner structures. It seems in 
particular, although information about the full stress 
tensor in the measurements of [l^, is lacking, that 
packings densified by vibrations or repeated shakes have 
a smaller coordination number than lubricated ones, in 
spite of a possibly larger density. Additional experimen- 
tal results with more detailed information on stresses and 
anisotropy of elastic moduli in packings assembled by dif- 
ferent procedures, as well as simulations of anisotropic 
packings, could enable more quantitative comparisons. 

Based on those results, several interesting perspectives 
should be pursued in the near future. 

On the theoretical side, it seems promising to study 
how granular packings, within and outside the elastic 
range, deform and destabilize, in more microscopic de- 
tail. Basically, packings with few contacts are closer to 
failure, and some of their anomalies in elastic proper- 
ties correlate with failure mechanisms. Amorphous sys- 
tems made of model atoms or molecules at zero temper- 
ature have been characterized, in this respect, in terms 
of an intrinsic scale [6ll . [6^ , and elementary plastic re- 
arrangement events, the spatial structure of which is 
similar to that of nonafRne elastic displacements, have 
been investigated (63j . Granular materials have friction, 
which requires more sophisticated stability analyses of 
given contact networks [H, [H, [H, H^, interact with 
much stiffer force laws, and exhibit characteristic dila- 
tancy properties and fabric evolutions under strain. It 
is worth investigating in greater detail the possible simi- 
larities and differences between their quasistatic rheology 
and the plasticity of amorphous materials. A character- 
istic length scale, which diverges in the isostatic limit, 
has also been invoked in relation with the singular elas- 
tic properties of frictionless packings [50|. It should be 
examined whether such a length plays a role in nonelastic 
deformation behaviors, and similar investigations should 
be carried out in systems with friction, which also exhibit 
complex, long-range correlated strain fields |62| . 



More practical issues which deserve investigations are 
how elastic moduli, which can be measured in equi- 
librated packings under varying stresses [13], can be 
used to infer useful information on their inner structure, 
which, in turn, can be exploited to predict their behavior 
under larger disturbances. As an example, the coordina- 
tion number of an isotropic packing, if it is large, will re- 
sult in a stiff response to deviator stress increments, char- 
acteristic of a stable contact network, and a faster mobi- 
lization of internal (macroscopic) friction as a function of 
strain (see Fig. [T7I and Ref. [5^). Numerical simulations 
of anisotropic stress and fabric states, of stress paths and 
large strains, and further numerically based correlations 
between elastic properties and stresses, strains and inner 
structures are of course necessary. Finally, the geometry 
of polydisperse and non-spherical particles should also be 
explored. Such packings might have large populations of 
rattlers, some more collective stable floppy modes than 
in the case of spheres f68l| , and different mechanical prop- 
erties 16911. 



APPENDIX A: TANGENTIAL ELASTICITY AND 
FRICTION 

We investigate here the effect on elastic moduli of the 
corrections to the contact law advocated in ^] that we 
use in our simulations, and we also discuss the possible 
effects of more sophisticated models, in which the par- 
tial mobilization of friction and the presence of a sliding 
region within the contact area [H, [7l| is taken into ac- 
count. 

Strictly speaking, those terms preclude the definition 
of a perfectly elastic response, which should be reversible 
and involve a uniquely defined stiffness matrix. 

The simplified law we adopt involves a tangential stiff- 
ness Kt depending on the normal deflection h, but in- 
dependent of the current mobilization of friction. This 
is the same approximation as used in [12] : the value 
of Kt is the correct one in the absence of elastic relative 
tangential displacement, when T = 0. 

The rescaling of Kt we chose to apply in situations of 
decreasing normal force TV, in order to avoid energetic 
inconsistencies, as explained in paper I [l[, means that 
the contact stiffness matrix, as defined in Section IIIBi 
should then be written differently. Its block correspond- 
ing to contact i,j is not given by IC^., as written in 

— 

but takes another form )C^. for a receding pair of grains. 
— *j 

Since Kt/ Kn is constant, and Kfq (x N^^^, one has, with 
and Ty/]]: 



as first and second basis vectors, 
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KT{h,_ 



(Al) 



The non-diagonal element of (jAip is smaller than 
^iK^ih)/?! (ifAf/lO for /i = 0.3), and its effects are likely 
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limited if friction is not too strongly mobilized, as should 
be the case under isotropic loads. Let us denote as AK 
the correction to the symmetric form of stiffness matrix 
^ j£ (|ggg Section III B|) . in which Eqn. ([6]) is applied 
to all contacts, due to this treatment of decreasing nor- 
mal forces. We solved the linear system of equations 
with loads corresponding to different global stress incre- 
ments, to first order in the perturbation AK: 

U ~ U(") + AU, 

where U'^"-' is the solution to the unperturbed problem, 
U(") = K-i . F°''*, and AU is the correction: 



AU = -K"^ • AK • U("\ 



(A2) 



In (jA2[) one should pay attention to the directional de- 
pendence of AK (according to whether h increases or de- 
creases). The first-order correction is therefore not lin- 
ear, but depends linearly on the amplitude of the load 
increment with a coefficient depending on its direction. 
We evaluated the resulting correction to the compliance 
in the cases of uniaxial {e.g., Aui > or Acri < 
while A(T2 — A(T3 — 0), isotropic (positive or negative 
value of AiJi = Ao'2 = Acra) and purely deviatoric {e.g., 
A(Ti — —Aa2 and Acts = 0) stress increments. Rela- 
tive corrections never exceeded 1%, the largest ones, as 
expected, being observed for an isotropic pressure reduc- 
tion (which tends to reduce normal contact forces). 

Our contact model also introduces an approximation, 
which we now discuss. As opposed to more sophisticated 
implementations of the contact law, as used by some au- 
thors [13, [3^ , our model does not keep track of the local 
slip distribution within the contact region. The maxi- 
mum effect of such slip is to reduce the tangential stiff- 
ness from Kt{N) to 



K't{N,T)=Kt{N){1 



N 



1/3 



(A3) 



in the "loading" direction [i.e., tending to increase 
||T||/A^), and for the tangential relative displacements 
along T. The possible influence on the simulated elastic 
properties of our overestimating the tangential stiffness 
of the contacts was assessed as follows. We computed the 
elastic moduli for the equilibrated configurations, keep- 
ing the same values of contacts forces, using formula [A3I 
for all contacts, in both tangential directions. (Such a 
calculation thus implicitly assumes that the equilibrium 
force distribution is not affected by the change in the con- 
tact law). This procedure exaggerates the effects of slip 
and gradual friction mobilization: formula IA3I gives the 
lowest possible value for Kt and only applies in specific 
loading histories, and for stress increments that tend to 

ITI I 

increase -^j^- We found that the relative corrections to 
computed elastic moduli evaluated with this procedure 
never exceed 3%. 

Consequently, it is a very good approximation to re- 
place the contact stiffness matrix IC by its diagonal form 



given by Eqn. ([6]), provided friction is not fully mobilized 
in any contact. 

If condition ||T|| = fiN is reached in contact then 
matrix JC has to be written as follows. With the same 

choice of basis vectors as for (jAip , 1C_^^ has a "loading" 

form JC^^ given by 



KN{h^j) 
fiKNih^j) 
Krih^j) 



(A4) 



and an "unloading" one equal to IC^_. or to ICv. , depend- 
ing on whether Sun is increasing or decreasing. If it is 
increasing, the loading form should be used if 



KTih,j)SuJ 



tiKN{h,j)du^^ > 



If Sun is decreasing, this condition becomes 



KT{h.,,)5uf^ ■ ^ 



KN{h^,)Suf^ > 0. 



Note that IC .^ is a non-symmetric singular matrix of rank 
2. As remarked in Section fll CI well-equilibrated config- 
urations prepared by molecular dynamics do not contain 
any contact where the condition ||T|| = /xiV is exactly 

reached. This justifies our using the diagonal form )C^. 

— u 

of the local stiffness matrix block pertaining to any con- 
tact. Full mobilization of friction occurs under small load 
increments, and the resulting deviation from elastic be- 
havior is investigated in Section IVIII 



APPENDIX B: GEOMETRIC STIFFNESS 
MATRIX AND TREATMENT OF DIVALENT 
BEADS 

The geometric term added to the change in intergran- 
ular forces entailed by small displacements and rotations 
was evaluated in paper I Here, for completeness, 
we write down its contribution to the geometric stiffness 
matrix K*-^-* . Those terms stem from the change in the 
direction of previous contact forces. They were carefully 
evaluated in general situations of particles of arbitrary 
shapes by Kuhn and Chang [g^I, and by Bagi This 
results in formulae of considerable complexity, involving 
"branch vectors" , joining contact points to particle cen- 
ters where moments are evaluated, as well as particle 
surface curvatures, which determine the small changes in 
normal directions caused by particle displacements and 
rotations. Those results assume much simpler forms in 
the case of spherical grains, for which both branch vec- 
tors and radii of curvatures are given by particle radii, 
and our results agree with such particular forms of the 
general expressions of Refs. [6^ [65|. 

In the matrix block Our results agree with the general 
expressions written down in those references, but radii of 
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sphere i is simply denoted as Ri below (in our numerical 
computations, all Ri values are equal to a/2). 

The 6x6 block K*-^-* of the geometric stiffness matrix 
is a sum over the contacts of grain i, 

each term being given by the expressions written in paper 
I [l|, appendix B] . Using a system of coordinates with riy 
and Ty setting the orientations of the two first axes, one 
has: 
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As to the non-diagonal block Kj^', it is obtained from 
~L^. on reversing the signs of the coefficients in the three 

first columns. Matrix K'^' is therefore clearly not sym- 
metric, which, in principle, forbids the definition of an 
elastic energy. However, each term involving Nij or 

in K*-^-* is negligible once compared to its counterpart in 

K^^\ where forces are replaced by terms of order Kjqa. 

Generally, K^^-* terms are therefore of relative order kT^ 

if compared to the corresponding ones in K*-^''. Thus 
the geometric stiffness matrix only plays a role for those 
directions of displacement vectors belonging to the null 
space of K'^'. This is important for frictionless grains, in 
which case such floppy modes of the constitutive matrix 
are necessarily unstable for spheres, but not so for, e.g., 
ellipsoids ^32, 68]. In the case of frictional spheres we 
did not obtain any floppy mode on the backbone, except 
for beads with two contacts. However, the correspond- 
ing free motion was shown to create no instability 1] . In 
order to work with a positive definite stiffness matrix, as 
rendered necessary by numerical techniques like Cholesky 
factorizations or conjugate gradient iterations, a stiffness 
term is added which associates an elastic energy with the 
free motion of divalent beads (the same trick is used to 
impede the global translations of all grains as one rigid 
body). Those particles consequently have no rotation in 
the direction defined by their two contact points, and the 
solution of the system of linear equations defined by (fT3|) 
is unique. 

With this precaution we can therefore safely neglect 
the geometric stiffness matrix in all cases studied in the 
present numerical work. 



APPENDIX C: ELASTICITY OF A 
PRESTRESSED SYSTEM 

We recall some basic properties which lead to a distinc- 
tion to be made, whenever the reference configuration of 
an elastic continuum is prestressed, between the elastic 
moduli and the coefficients of the linear relation between 
the Cauchy stress tensor and the strains. We specify 
here which type of elastic coefficients we compute in sim- 
ulations. Then, we also point out that some small cor- 
rections -negligible in practice - should in principle be 
applied to the computed moduli, because of preexisting 
stresses. 

Let us consider a uniform displacement gradient Vu 
within an elastic continuous medium (derivatives with 
respect to coordinates on the reference, undisturbed con- 
figuration). In linear elasticity, the free energy density, 
A/VLq, evaluated in a reference configuration (with vol- 
ume r2o)j is a quadratic function of the Green-Lagrange 
strain tensor e, which expresses material deformation [1[ . 
The first-order term is written with the Piola-Kirchoff 
stress tensor tt^ [t^ in the reference configuration: 



(CI) 



To second order, the free energy associated with small 
strains involves the tensor of elastic constants C: 



A = no 



lo ■ i + 2= ■ = • = 



Elastic moduli thus appear in the increment of tt: 



TT — TTj^ = C : e. 



(C2) 



(C3) 



The Voigt symmetry Cap-yS = G^sap is satisfied by 
the coefficients of this linear law because they are sec- 
ond order derivatives. Let L denote the diagonal ma- 
trix containing the cell dimensions along the three axes, 
and its value in the reference configuration, corre- 
sponding volumes being f7 and fio- With the notation 
£ = 1 -I- Vu = ^ • one has n/^o = detg for the 
dilation and g_ is related to tt by 

^= (det£)£^^ -^^ "^^-^ (C4) 

(This relation between the Piola-Kirchoff stress tensor 
and the Cauchy stress tensor can be found in continuum 
mechanics textbooks [t^ and was recalled in paper I). 

Let us now use (jCSP and (|C4p to write the increment 
of the Cauchy stress tensor to first order in displacement 
gradient Vu (needless to distinguish spatial derivatives 
in the reference or the deformed configuration at this 
stage, as this would introduce second order corrections). 
Defining the linearized strain tensor e as 



1 



( Vu 



'Vu), 
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one obtains: 

g^Eo~ lo ~ = ' lo ~ lo ' ^= + = • i- 

Therefore, a is not necessarily a function of the symmet- 
ric part of Vu only. A rigid rotation (for which Vu is 
antisymmetric) might produce a Cauchy stress increment 
if TT^ and Vu do not commute. Likewise, the coefficients 
expressing the linear dependence of a on Vu do not al- 
ways satisfy the Voigt symmetry, and hence one cannot 
regard a constant a as deriving from a potential energy 
of external loading. Both conditions, symmetry and de- 
pendence on e only, are however restored if one restricts 
to symmetric displacement gradients, or if n and e share 
common principal directions. This is always the case 
with our choice of boundary conditions, or in general if 
TT^ = PI is an isotropic tensor. In this latter case, (jCSP 
relates g_ to e, assuming isotropy of the material, with a 
tensor of "apparent" elastic moduli B, that has the same 

symmetries as C. C, in isotropic systems, can be writ- 
ten with a bulk modulus B and a shear modulus G. On 
relating a to e, the apparent moduli (as measured in an 
experiment) are B + P/3 and G — P. 

It should be specified that our procedure to compute 
elastic moduli in isotropic sphere packings is based on 
a formula for the Cauchy stress tensor. Therefore, the 
resulting moduli are the elements of matrix B, rather 

than C. 

As a consequence of the stress and forces that pre- 
exist in the initial configuration before elastic response 
is probed, our results should also in principle be slightly 
modified. As we use it, the equilibrium equation for stress 
components actually gives the increment of the product 
QcTq, , from which the contribution Aficr^ , due to volume 
change Afl = — figtre should be subtracted before divid- 
ing by if the Cauchy stress variation is to be obtained. 
As a consequence of this correction, we should add P/3 to 
the value of B obtained with our calculation procedure. 
This is a very small effect, which we have been neglecting 
(moduli are in MPa for stresses in kPa, see Fig. [1]). 

APPENDIX D: VOIGT AND REUSS BOUNDS 
FOR ELASTIC MODULI IN A SPHERE 
PACKING 

Within the approximation that the stiffness matrix 
does not depend on the direction of the stress (or strain) 
increment, and is symmetric (see Section III CI and Ap- 
pendix [X| , which fortunately proves accurate (see also 
Section fVIip . the elastic regime can be defined, and the 
variational properties (fT5|) and (fT6|l can be used. Vari- 
ational properties leading to bounds for moduli are sel- 
dom invoked in the context of granular materials. Our 
purpose here is to recall how these useful properties are 
established and interrelated, and how they can be ex- 
ploited. 



Let us first briefly establish the minimization property 
for contact force increments, which is less familiar than 
the simple minimization of elastic energy (|15p . We con- 
sider the solution Af * to the problem of minimizing (jl6p 
among contact force increment vectors that balance the 
applied load increment, i.e. such that 

• Af AF°'^* (Dl) 

This solution is characterized by the existence of a vector 
X of Lagrange multipliers such that 

• Af* £ • X, 

and (|Dip thus entails 

'^£-^-£-x=K-x = AF*^=^', 

This means that x is the displacement vector solution to 
the elastic problem, and that Af * = £ • G • x is indeed 
the corresponding contact force increment vector. 

We now derive the explicit formulae for Voigt and 
Reuss bounds. 

A first step is to exploit the isotropy of the medium, 
which enables inequalities to be written separately for 
bulk and shear moduli. In our stress-controlled approach, 
Act is imposed, minimum values in p5|) and (|16|) are 

= -^A^:£-^ : A^ 

The values obtained with trial solutions for displacements 
or contact force increments can then be regarded as es- 
timates of those quadratic expressions in Act, and hence 
provide estimates of the corresponding compliance ma- 
trix C~^. The meaning of C^^, in a finite sample, is 

specific to the choice of particular boundary conditions. 
In the large sample limit, it is assumed to satisfy the 
symmetry properties of the medium, which is statisti- 
cally isotropic in the numerical studies reported here. 
Moreover, it is also expected to approach the macro- 
scopic compliance matrix, whatever the particular choice 
of boundary conditions. If, as in our numerical study, we 
restrict Act to a diagonal form, and hence regard it as a 
vector with three coordinates CTq, a = 1,2,3, is a 
matrix S of the form 





S12 


Si2 




Sii 


Si2 


S\2 


S12 


Sii 



with, due to isotropy, 

^'''9B^3G 
9B 6G' 
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When Aa is an isotropic pressure increment AP, one 
has 



»o(AP)^ 
3B 



whence an upper bound to B with an estimate of W^, 
and a lower bound with an estimate of W2 ■ 
When Act is of the form (g, —q, 0), then 



2G 



(D3) 



whence an upper bound to G with an estimate of W*, 
and a lower bound with an estimate of W2 ■ 

The Voigt bounds on B and G are based on trial dis- 
placements defined as (see Eqn. [S]) 



U = ((0, 0)i<j<„, (eQ)l<a<3) 



(D4) 



The choice of e should then minimize psp . restricted to 
displacements of this particular form, i.e., 



no 



in which the 3-vector notations for strains and stress in- 
crements as defined above are adopted, and 3x3 matrix 
L is a sum over contacts, to be evaluated below. The 
best choice for e is 



e* = • CT 



The minimum value of Wi then yields 

gVoigt _ j^-l 



(D5) 



Let us now write down matrix L in a sphere packing with 
our boundary conditions. We introduce the notations 



T 



L. 



for each contact i,j between spheres of radii Ri and Rj, 
and neglect hij in comparison with the radii, as we have 
been doing throughout this article. Then one has, for 
each pair of indices a, (3, 1 < a,(3 < 3 (no sum over 
repeated indices) 



-y 

i<3 



(D6) 

In (|D6p . 6ap is the Kronecker symbol and index pairs 
i < j run over the list of force-carrying contacts between 
grains labelled i and j. 

In general, it might be favorable to allow for a com- 
mon spin, i.e., A9i — Co, for all particles in the choice 
trial displacements (jD4[) . The optimal choice uj* only 
vanishes when the stress tensor and the fabric tensor de- 
fined as the average of n^j ® n^j over contacts, weighted 
by LJ, , share the same eigenvectors. This conclusion, wich 



was reached before by Jenkins and La Ragione [73|, and 
independently by Gay and da Silveira [7^], on directly 
estimating the stress increments corresponding to a pre- 
scribed strain, is retrieved here as an illustration of the 
variational approach. 

Returning now to the case of isotropic packings of 
monodisperse spherical beads of diameter a, the fabric 
tensor is isotropic, which ensures lo* . Thus matrix L is 
the Voigt estimate of C. To compute its terms in the 
large system limit for isotropic packings, we transform 
the sums in (jD6|) as we did for the evaluation of average 
stiffness by (IT51) . Exploiting the independence between 
stiffness fluctuations and contact orientations, as well as 
relations 

{ni) = ^ and (n^n^) = ^ 
for isotropically distributed unit vectors, one gets 

C-^* = ^ f V'' ^±^Z(1/3)PV3 

2 \ TT / 15 



^Voigt_ 34/3 ( z^Ey \-aT 



Z(l/3)pi/3^ 



from which expression of pV°'g' and G^"'*^! jg readily 
derived, since B = (Cii-h2Ci2)/3 and G = (Gii-Gi2)/2. 

Finally, to establish the Reuss lower bound for B, we 
choose an isotropic stress increment Act = (AP, AP, AP) 
and evaluate W2 for a trial set of contact force increments 
chosen, in any sample in equilibrium under pressure P, 
as 



Af,. 



AP^ 

P 



(D7) 



in contact i,j, initially carrying force f^ . Such force in- 
crements balance the load increase AP by linearity of 
equilibrium relation (|Dip . The resulting value of W2, 



W2 = 



is quadratic in AP, and yields B^°^^^ , as written in 
on identifying the optimal 14^2 value with the macroscopic 
energy: 



w 



once the sum is transformed by the same procedures as 
in the evaluation of the average stiffness (Kn) in (|19|). 
using the definition of Z{5/3) in 

No such trial vector of contact force increments as (|D7|) 
is readily available when the applied stress increment 
is not proportional to the preexisting stress, which is 
isotropic in the present study. In general, in anisotropic 
stress states, a similar Reuss-type approach is expected 
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to provide a lower bound estimate for a certain combi- 
nation of elastic moduli, which expresses the response to 
proportional load increases. 

Finally, let us recall that similar minimization prop- 
erties as (fTS)) and (fT6)l hold in a strain-controlled ap- 
proach. If strains e are imposed, then displacement vec- 
tor U, which is constrained to correspond to e (this sets 
the values of its three last coordinates with our choice of 
boundary conditions) should minimize: 

5i(U) = iu-K-u 



while the contact force increments, vector Af, should 
equilibrate each grain and minimize 

£'2(Af) = ^Af • • Af - ilAg : g, 

in which the stress increment Aa is directly given by Af , 
just like stress components relate to contact forces in ([T]). 

It is easy to check that the strain-controlled approach 
yields exactly the same Voigt and Rcuss bounds as the 
stress-controlled one. 
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